Non-monotonic residual entropy in diluted spin ice: a comparison between Monte Carlo 
simulations of diluted dipolar spin ice models and experimental results 
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Spin ice materials, such as Dy2Ti207 and Ho2Ti207, have been the subject of much interest for over the past 
fifteen years. Their low temperature strongly correlated state can be mapped onto the proton disordered state 
of common water ice and, consequently, spin ices display the same low temperature residual Pauling entropy 
as water ice. Interestingly, it was found in a previous study [X. Ke et. al. Phys. Rev. Lett. 99, 137203 
(2007)] that, upon dilution of the magnetic rare-earth ions (Dy 3+ and Ho 3+ ) by non-magnetic Yttrium (Y 3+ ) 
ions, the residual entropy depends non-monotonically on the concentration of Y 3+ ions. In the present work, 
we report results from Monte Carlo simulations of site-diluted microscopic dipolar spin ice models (DSIM) that 
account quantitatively for the experimental specific heat measurements, and thus also for the residual entropy, 
as a function of dilution, for both Dy2-a: Y^T^C^ and Ho2- 2 ,Y :c Ti207. The main features of the dilution 
physics displayed by the magnetic specific heat data are quantitatively captured by the diluted DSIM up to, 
and including, 85% of the magnetic ions diluted (x = 1.7). The previously reported departures in the residual 
entropy between Dy2- a: Y I1 Ti207 versus H02-1 Y^^Cv, as well as with a site-dilution variant of Pauling's 
approximation, are thus rationalized through the site-diluted DSIM. For 90% (x = 1.8) and 95% (x = 1.9) of 
the magnetic ions diluted, we find a significant discrepancy between the experimental and Monte Carlo specific 
heat results. We discuss some possible reasons for this disagreement. 

PACS numbers: 75.10.Hk, 05.50.+q, 75.40.Mg, 75.50.Lk 



I. INTRODUCTION 

The theoretical and experimental study of geometrically 
frustrated magnets 1 - constitutes a very active research area in 
contemporary condensed matter physics. In these systems, the 
predominant interactions compete with each other, inhibiting 
the development of long-range magnetic order down to very 
low, if not zero, temperature.- The temperature regime where 
strong magnetic correlations exist, but long range order is ab- 
sent, is commonly referred to as spin liquid 4 or cooperative 
paramagnetic stated 

One particularly topical example of geometrically frus- 
trated magnets is spin ice materials £2tlr— These are realized 
by the canonical compounds Dy2Ti207 and H02TI2O7, as 
well as by the less extensively studied Dy2Sn20y — and 
Ho2Sn20yjii More recently, high pressure chemical synthe- 
sis has allowed one to make the Dy2Ge207 and Ho2Ge207 
compounds, and thermodynamic measurements have shown 
these materials to be an interesting new class of spin ice sys- 
temsJ^ii 3 - In that context, it is interesting to note that the 
CdEr2Se4, in which Er 3+ is unusually described by an Ising 
spin, as also been shown to be a spin ice. 14 » In all of these ma- 
terials, the magnetic rare-earth Dy 3+ , Ho 3+ and Er 3+ ions sit 
on the vertices of a pyrochlore lattice of comer-sharing tetra- 
hedra; the Ti 4+ , Sn 4+ and Ge 4+ ions are non-magnetic. Be- 
cause of the very large single-ion anisotropy at play in these 



systems, the moments can be described below a temperature 
T ~ 50 K as classical Ising spins pointing along the local 
[111] direction at their respective pyrochlore lattice sites. 8,15,16 
Below a typical temperature of order 1 K, the magnetic state 
of (Dy,Ho)2(Ti,Sn,Ge)207 can be mapped onto the proton dis- 
ordered state of common water ice^ hence the name spin 
ice£ In this low temperature spin ice state, the magnetic mo- 
ments are highly correlated locally and obey the so-called "ice 
rules": two spins point in and two spins point out of each tetra- 
hedron of the pyrochlore lattice, but without displaying long 
range order^ The spin ice state can thus be viewed as a co- 
operative paramagnet, 6 * or a classical spin liquid to adopt a 
more modern terminology* 4 - The label "classical spin liquid" 
stems from the very strong Ising nature of the lowest-energy 
crystal-field doublet for both Dy 3+ and Ho 3+ which results 
in a dramatically quenched level of quantum spin dynamics, 16 
At the same time, the high energy barrier to single spin flips 
causes the relaxation dynamics to become very slow in these 
materials below T ~ 1 K. Consequently, spin ices should be 
viewed as extremely sluggish classical spin liquids J£ 

For water ice, extensive calorimetric studies had been car- 
ried out long befor e - 19 ' 20 its magnetic counterparts were dis- 
covered^ The nature of the proton disorder in ice was de- 
scribed by Linus Pauling who estimated the residual en- 
tropy to be Sp = R/2m(3/2) per mole of protons (R 
is the molar gas constant), 1 - matching closely with experi- 
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mentSjii2P. The same residual entropy is found in spin ice ma- 
terials j 9 ' 10 ' 13,14 ' 21,22 providing direct thermodynamic evidence 
for the ice rules being at work. In spin ices, the crossover 
from the paramagnetic phase to the macroscopically degen- 
erate spin ice state with its Pauling residual entropy is sig- 
nalled by a broad specific heat peak at around T p ~ 1.2 K 
for Dy 2 Ti 2 OT^3 and T p - 1.9 K for Ho 2 Ti 2 7 ^ There is 
no thermodynamic phase transition between the high temper- 
ature paramagnetic state and the low temperature spin ice state 
as evidenced by the absence of sharp thermodynamic features 
at T p . Theoretical studies have shown that long range mag- 
netostatic dipole-dipole interactions are responsible for the 
finite entropy spin ice state in (Ho,Dy) 2 (Ti,Sn,Ge) 2 07 com- 
pounds.— _ — Yet, it is generally theoretically accepted that the 
same dipolar interactions should give rise to long-range order 
at a critical temperature T c <c T p if true thermal equilibrium 
could be maintained down to sufficiently low temperature.^^ 
Indeed, Monte Carlo simulations that employ loop moves to 
generate non-local spin flips, which allow the system to re- 
main in thermal equilibrium without violating the ice rules, do 
find a transition to long-range order at low-temperature^ 2 ^ 
To this date, however, no experiment has found a transition to 
long-range order in spin ice materials^ presumably because 
of a dynamical arrest in spin flips 3 -^ and the associated relax- 
ation times growing exponentially fast below a temperature of 
about 1 K. 

Considering the broader context of ice-like systems devel- 
oping extremely slow dynamics at sufficiently low tempera- 
tures, one notes that as water ice is doped with alkali hydrox- 
ides, such as KOH or RbOH, a sharp first order transition to 
long range order occurs at a temperature near 72 K. At that 
transition, a large portion of the residual Pauling entropy is 
released through the latent heati 3 ^ 3 . These experiments sug- 
gest that the proton-disordered ice state is somewhat fragile 
against impurities and that the frustrated disordered ice state 
with residual entropy can be eliminated through the influence 
of impurities and/or random disorder. Yet, despite much theo- 
retical work, it remains unclear what is the precise mechanism 
via which alkali hydroxides in the water ice system promotes 
the development of long range order*2i 

Inspired by the impurity-driven long-range order observed 
in water ice* 3 ^ 3 . it is interesting to ask whether the mag- 
netic spin ice analogue could also display interesting behavior 
when subject to the addition of random impurities. For ex- 
ample, perhaps a slight dilution of the magnetic Dy 3+ and 
Ho 3+ ions could lower the kinematic barriers for spin flips, 
thus accelerating the spin dynamics, and help promote a tran- 
sition to long-range order without significantly affecting the 
broken discrete symmetry long-range ordered ground state 
of dipolar spin \czM^ In that context, we note that mag- 
netic site-dilution in spin ices can be realized rather straight- 
forwardly in the Dy 2 _ a; Y 2 ;Ti 2 07 and Ho 2 _ a; Y z: Ti 2 07 com- 
pounds, which form a solid solution over the whole x G [0, 2] 
range, and where the magnetic Dy 3+ and Ho 3+ ions are re- 
placed by non-magnetic Y 3+ ions™ 3 - The close ionic radius 
of Y 3+ with that of Dy 3+ and Ho 3+ allows for a substitution 
that causes negligible local lattice deformation and strain. Di- 
lution of Dy 3+ /Ho 3+ by Y 3+ can thus be viewed as a mere 



replacement of the Dy 3+ /Ho 3+ magnetic species by a mag- 
netically inert substitute. A recent neutron scattering experi- 
ment shows no sign of long-range ordering in Ho2- x Yj;Ti 2 07 
down to 30 mK for x = 0.3 and x = 1.0^ On the other 
hand, specific heat measurements have found that the low- 
temperature residual entropy, S Tes , of diluted Dy 2 _ x Y x Ti 2 07 
and Ho 2 _ 2; Y 2: Ti 2 07 spin ices display a non-monotonic de- 
pendence on the level of dilution^ A calculation generaliz- 
ing Pauling's argument (gPa) to the case of site dilution of 
a nearest-neighbor spin ice model^ was able to qualitatively 
account for such a non-monotonic behavior. 23 However, the 
apparent systematic departures between the gPa and the ex- 
periment results as well as the differences between Dy- and 
Ho- based materials (see Fig. [TJ have so far remained un- 
addressed. It was suggested in the original work.— that the 
residual entropy may be material-dependent and have a more 
drastic non-monotonic dependence on levels of dilution than 
the analytic generalized Pauling argument (gPa) does. The 
reason for these differences might be caused, for example, 
by the extra complexities of the long-range dipolar interac- 
tions compared with the nearest-neighbor model. In this pa- 
per, we address and rationalize quantitatively the origin of 
the difference in residual low temperature entropy between 
Dy 2 -a;Y x Ti 2 07 and Ho 2 _ x Y 2; Ti 2 07 as well as with the gPa 
illustrated in Fig.Q] 
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FIG. 1: (Color online) Adapted from X. Ke et. al. Phys. Rev. Lett. 
99 137203 (2007). Experimental residual entropy as a function of 
dilution level x (see Ref. I23T0 . Dy denotes Dy2- a: Y a: Ti207, Ho 
denotes Ho2- a; Y :E Ti207 and Gen. Pauling denotes the generalized 
Pauling approximation (gPa) presented in Ref. I23I1 . As noted by Ke 
et. al, there is an obvious systematic departure between the three 
curves, except for the undiluted compounds (x = 0). 



Moving away from the specific context of disorder and im- 
purities in ice-like (water or spin) systems, one notes that 
the problem of quenched random disorder in highly frustrated 
magnetic systems is one of long-standing interest, going back 
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to the seminal work of Villain,* In more recent years, the ef- 
fects of disorder on the thermodynamic properties of highly 
frustrated magnetic systems, in large part motivated by studies 
on kagome materials such as SrCr 2 .Gai2_ a: Oi9 (SCGO)^— 
and ZnCu3(OH)gCl2 (Herbertsmithite)^^d£ has been a topic 
of much interest. Research efforts in this area have been es- 
pecially motivated by the necessity to understand whether the 
observed experimental behavior in SCGO and Herbersmithite 
is intrinsic to the hypothetical disorder-free material or is, in- 
stead, (largely) driven by impurity effects. The problem of di- 
lution in quantum triangular antiferromagnets is also one with 
very interesting and rich physics4i^ 

Spin ice is at the present time one of the best understood 
highly frustrated magnetic systems, both from a microscopic 
model perspective 4 - 3 - as well as from a field theory one.^ 4 -^ 
Spin ices would thus appear to be an ideal system to investi- 
gate quantitatively the effects of random disorder in a highly 
frustrated magnetic setting.^ 4 ^ This is precisely the broader 
goal of this paper: to perform such a quantitative compari- 
son between theoretical modeling and experimental measure- 
ments in a specific class of disordered highly frustrated mag- 
netic materials. As a first agenda in this program, we con- 
sider the aforementioned problem of diamagnetic site-dilution 
dependence of the residual entropy in Dy2- a; Y 2 ;Ti207 and 
Ho2- 2; Y 2 ;Ti207 spin ice materials. 

It turns out that there is a growing interest in the problem of 
disorder in magnetic pyrochlore oxides. For example, direct 48 
and indirect 49 evidence has recently been put forward that, 
in image furnace grown single crystals, there is a small level 
(0(1%)) of substitution of the Ti 4+ transition metal ions by 
trivalent rare-earth ions — a phenomenon referred to as "stuff- 
ing". Other examples include the mixing of different types 
of ions on the rare-earth site^ or different non-magnetic ions 
at the B site.^i^ Thus, in comparison with these various dis- 
order settings, which would all generate random bonds, the 
problem of site-dilution may be expected to be simpler, and a 
necessary first step in our goal of understanding the effects of 
random disorder in magnetic pyrochlores oxides.— 

In order to investigate the microscopic origin of the relative 
departure of the three curves in Fig. Q] we performed Monte 
Carlo simulations of a diluted variant of the pertinent micro- 
scopic dipolar spin ice model of Hc^^Ot^ 4 - and Dy2Ti2074^ 
A direct comparison of the temperature-dependent magnetic 
specific heat, C m (T), for various dilution levels, x, between 
simulations and experiments is made in order to validate a 
simple site-diluted version of the otherwise pure (dilution- 
free) microscopic models. Through the simulation data, we 
obtain an accurate C m (T), which provides for a precise deter- 
mination of the residual entropy, down to the lowest temper- 
ature To <~ 0.4 K considered in experiments^ 3 The simula- 
tion results confirm the previous speculation^ 3 , that the depar- 
ture of the material-dependent residual entropy from the gen- 
eralized Pauling argument (gPa) occurs because of material- 
specific details of the interactions. That said, our conclu- 
sion regarding the difference in residual entropy S Ies be- 
tween Ho2Ti2C<7 and Dy2Ti2C<7 is different from the one in 
Ref. ll23"ll . namely, we find Ho2-^Y a; Ti207 to have a smaller 
S les (x) than Ho2- x Y x Ti207 does. 



The rest of the paper is organized as follows: In Section II, 
we discuss the details of the experimental methods; In Section 
III, we present our microscopic models and the Monte Carlo 
simulation methods; in Section IV, we present and discuss the 
results of the Monte Carlo simulations and address the previ- 
ously reported^ 3 , material-dependent residual entropies along 
with their departure from the gPa predictions; Section V con- 
cludes the paper. 



II. EXPERIMENTAL METHODS AND RESULTS 

Specific heat measurements were performed on Y-diluted 
spin ice materials, Dy2- x Y x Ti207 and I-k^-xY^^C^, us- 
ing a Quantum Design Physical Property Measurement Sys- 
tem (PPMS) cryostat with the He3 option via a standard semi- 
adiabatic heat pulse technique. The Dy-based samples were 
thoroughly mixed with Ag and pressed into pellets to facili- 
tate thermal equilibration. The scaled Ag specific heat, mea- 
sured separately, was subtracted from the total specific heat. 
The phonon contribution was extracted by fitting the data with 
the Debye formula in the temperature range T <E [10, 20] 
K, and subtracted from the total specific heat to obtain the 
magnetic specific heat contribution, C m (T). Ho-based sam- 
ples were pressed directly into pellets and the magnetic spe- 
cific heat was obtained after subtracting both the phonon and 
the large Ho nuclear Schottky anomaly contribution.^^ The 
data, C m (T)/T, integrated from T (x) = 0.4 ± 0.1 K, de- 
pending on the lowest temperature To (a;) accessed for a given 
concentration x, up to a ('high') temperature T, was used 
to determine the residual low-temperature entropy, S les (To). 
The previously reported^ residual entropy is reproduced here 
in Fig.[T]for convenience. As discussed in the Introduction, the 
residual entropy plotted in Fig. 1 varies non-monotonically as 
a function of the Y concentration for both the Dy2_a; Y^^C^ 
and the Ho2- x Y x Ti207 series, being qualitatively captured 
by a generalization of Pauling approximation's (gPa) that is 
represented by the dashed curve.— 



III. MICROSCOPIC MODELS AND MONTE CARLO 
SIMULATIONS 

A. Microscopic Models of Spin Ices 

In spin ices, the magnetic moments reside on a pyrochlore 
lattice, which consists of a face-centered cubic lattice of 
corner-sharing tetrahedra primitive units.^ Due to the large 
energy scale (~ 300 K) of the crystal field splitting between 
the ground state doublet and the lowest-energy excited dou- 
blet that exist in Dy 2 Ti 2 07 and rk^^Or , 5 ' 15 ' 16 the states that 
form the ground doublet of the Dy 3+ and Ho 3+ ions can safely 
be assumed to be the only thermodynamically relevant states 
below a temperature T < 50 K. 

As suggested originally^ the minimal model that describes 
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the geometrical frustration in spin ices is 

%NNSIM = Jefi ^ CjOj (1) 

where J e g > is the effective antiferromagnetic interaction 
between the er's Ising variables. This model possesses a Paul- 
ing residual entropy, Sp,— and displays at zero temperature 
an ice-rule obeying ground state characterized by dipolar- like 
spin-spin correlations that emerge from the "two-in"/"two- 
out" ice rule constraints^ 

On the other hand, in the real spin ice materials, the Dy 3+ 
and Ho 3+ ions carry a large magnetic moment (~ 10 //b) 
and the long range dipolar interactions cannot be ignored^2^ 
Given the symmetry of the crystal field ground state^i^ the 
magnetic moments can be well described by vector spins con- 
strained by the single-ion anisotropy to point strictly parallel 
or antiparallel to their respective local [111] direction (i.e., 
along the line from the corners to the centre of each tetrahe- 
dron) i 7,8 ' 15 Taking the dipolar interaction and the essentially 
infinite local Ising anisotropy into consideration, the dipolar 
spin ice model (DSIM) is defined by the Hamiltonian: 

r 3 

%DSIM = ^ SiSj < 2^ ^ v ^nf,r„ %i ' z j 

+ D(ri/r i3 f [zi ■ ij - 3 {z t ■ hj){zj ■ fij)} |(2) 

where Oi = ±1 are the Ising spin variables. The first term 
describes the Ising exchange interaction while the second term 
is the long-range magnetic dipole-dipole interaction. Here, 
v = 1, 2 or 3 refers to first, second or third nearest neighbors 
respectively, where J„ is the exchange coupling and r u is the 
distance between them. There are two types of third nearest 
neighbor interactions which we do not differentiate 4 s - Zi is the 
local [111] direction of the Ising axis and D is the strength of 
the dipolar interactions at nearest-neighbor distance. 

Using the most up-to-date values for J u and D that we 
are aware of, we have with our sign convention of the J v 's 
(J„ > is antiferromagnetic; J v < is ferromagnetic): 
Ji = 3.41 K, J 2 = -0.14 K, J 3 = 0.025 K and D = 1.32 
K for Dy^Y^O^ and J x = 1.56 K and D = 1.41 K 
for Ho2_ x Y x Ti207i21 Unfortunately, because of the complex- 
ity introduced by the large hyperfine coupling interactions in 
Ho-based materials, much less systematic calorimetric mea- 
surements, which provide many of the constraints to deter- 
mine Ji and J2r^ have been carried out on Ho2Ti2C>7 com- 
pared to Dy 2^07. Consequently, the J2 and J3 values for 
Ho2Ti207 have not yet been determined^! and we therefore 
set J2 = J3 = for this compound. As we shall see be- 
low, it turns out that this (J2 = J3 = 0) model describes 
well the magnetic specific heat of F^-zY^^C^ for the 
x = 0, 0.4, 0.8 and 1.2 values considered in this work. 

For the diluted samples, we assume that the non-magnetic 
diluting Y 3+ ions are introduced randomly while all other pa- 
rameters of the material, and therefore those of the model in 
Eq. (O, are assumed to be unchanged. This means that, un- 
til more accurate microscopic ab-initio modeling of the effect 



of diamagnetic site-dilution in spin ice compounds becomes 
available, we ignore local lattice strain effects that may re- 
sult from the substitution of Dy 3+ or Ho 3+ by Y 3+ . In prac- 
tice we thus ignore any changes that may occur in the J u ex- 
change couplings and the rare-earth ion magnetic moment /i 
that would result from variation of the single-ion crystal field 
ground state wavef unctions. This would seem a reasonable 
first approximation given the close ionic radius of Y 3+ with 
Dy 3+ and Ho 3+ . We note in passing that such an approxima- 
tion has recently been shown to describe quantitatively quite 
well the variation of the critical ferromagnetic temperature in 
Ho 3+ substituted by Y 3+ in LiHoi_ 2: Y a; F4 all the way tori^ 
and perhaps even including, the dipolar spin glass regime^S 
In practice, the microscopic J^'s and D in Eq. © are kept 
to their pure Dy2Ti2C<7 and Ho2Ti207 values while the Ising 
variables are redefined as Oi — > e^Ui, with e.; = if site i is 
occupied by non-magnetic Y 3+ ion or e, = 1 if occupied by 
a magnetic rare-earth ion. Thus, for [Dy,Ho] a; Y2_ a ;Ti207 the 
site-random probability distribution of e,;, -P(ei), is given by 
P( £i ) = (x/2)8(e t ) + (1 - x/2)5(ei - 1), where S(u) is the 
Dirac delta function. 



B. Monte Carlo Methods 

We carried out Monte Carlo simulations for the above 
model for Dy2_ K Y a; Ti207 and Fk^-xYz^CV at various 
Y 3+ concentrations x. We used a conventional cubic unit cell 
containing 16 spins, with the system of linear size L having 
16L 3 spins. Dilution is treated by randomly taking spins out 
of the system, and a disorder average over 50 different ran- 
dom dilution configurations was performed for each dilution 
level x. Periodic boundary conditions are used, and we imple- 
ment the infinite dipole interactions using the Ewald summa- 
tion technique^ Most of the data production was done with 
L = 4 while, for higher dilutions (x > 1.5), we used L = 5 to 
have a reasonably large number of spins remaining in the sys- 
tem. For most of the results presented below, very little system 
size dependence for the magnetic specific heat, C m (T), data 
was observed. 

A conventional single spin-flip Metropolis algorithm was 
employed for the Monte Carlo simulation. In addition, we 
used a non-local "closed-loop" updat o 28,29 as well as a new 
"open-loop" update that we now explain. The open-loop up- 
date is a modified version of the closed-loop update with the 
following amendments. In a diluted system, a fraction of the 
elementary tetrahedral units will have one or three sites occu- 
pied by a spin. Such "± tetrahedra" will have the sum of the 
Ising (Ji variables over the occupied sites equal to ±1 or ±3. 
At low temperatures, almost all such tetrahedra become con- 
strained to ±1, since these states are energetically lower than 
the ±3 ones. 

The open-loop update algorithm searches for an end-to-end 
chain of spins connecting two of these tetrahedra with oppo- 
site sums of the Ising variables. An open-loop update flips 
all the spins along the chain when accepted. Energetically, 
the nearest-neighbor part of the Hdsim is unchanged in such 
an open-loop Monte Carlo update. We use the term open- 
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loop update to stress the similarity of the algorithm to the 
original closed-loop update*2&2iL but with the chains of the up- 
dated spins ending at two "±1 tetrahedra". In order to further 
facilitate the equilibrium of the system, we found it neces- 
sary to also employ the parallel tempering technique which 
is commonly used in the study of spin glass models.— At 
least 200,000 Monte Carlo update steps are used with each 
single-spin-flip update sweep followed by the two types of 
loop moves update as well as by a parallel tempering replica 
exchange sweep. 5 - Another 200,000 such steps are used for 
data production. The magnetic specific heat was determined 
by performing a disorder average of the energy fluctuations: 

C (T) [{E2) {E)2] (3) 
CmU } " k B T 2 (3) 

where (...) and [...] are thermal and disorder averages, respec- 
tively. 



IV. RESULTS AND DISCUSSIONS 

We plot in Fig. [2] the magnetic specific heat versus tem- 
perature, C m (T), obtained from Monte Carlo simulations 
of Eq. (O (solid lines) for various levels of dilution in 
comparison with experimental data (open black circles for 
Dy2- X Ya;Ti207, open red squares for Hc^-aY^^C^). 

The agreement between our Monte Carlo simulation and 
the previous experiment^ is strikingly good for most dilution 
levels (up to and including x = 1.7 for Dy2- ;c Y 2 ;Ti207) and 
over a rather wide temperature range T ~ [0.4 K — 5 K]. This 
is particularly noteworthy given that there is no adjustment of 
the microscopic parameters of the dipolar spin ice Hamilto- 
nian of Eq. (O, except for the dilution of spins in the system. 
From these results, we can immediately conclude that a sim- 
ple site-diluted version of the DSIM of Eq. (fJI does capture 
the dilution physics of both materials at a quantitative level. 
This constitutes the main conclusion of this paper. 

Close inspection of Fig. |2]shows that there is a discrepancy 
in C m (T) between simulation and experimental results for 
T > 5 K. Also, the simulation results show a rise of C m (T) as 
T decreases below a temperature of approximately 0.4 K and 
0.6 K for Dy 2 _ 2: Y :! .Ti207 and Fk^-^Y^^Cv, respectively, 
while this behavior is barely noticeable in the experimental 
results. We address these two points in further detail in Sub- 
section IIV Al mostly at the phenomenological level, postpon- 
ing the discussion of the physical implications of these results 
for the determination of the residual entropy in the following 
subsection. In Subsection IIV Bl we present the low tempera- 
ture limit {To) dependence of the residual entropy, S ICS (To), 
as a function of dilution level, x. We comment in Subsection 
IIV CI on the failure of our Monte Carlo simulations to repro- 
duce the experimental results for x = 1.8 and x = 1.9. 



A. High and Low Temperature Regimes 

1. High temperature regime 

In the "high-temperature regime", typically above 4 K ~ 
5 K, we observe that our simulation results for C m (T) de- 
part from the experimental results. Such discrepancies need 
clarification since (i) a demonstration of the validity of the 
microscopic models considered depends on achieving a good 
degree of agreement between experimental and Monte Carlo 
C m (T) curves and since, (ii) as we shall see when discussing 
the residual entropy in the next subsection, C m (T) for T > 5 
K contributes up to about 10% of the full Rln(2) magnetic 
entropy. 

From a high-temperature expansion perspective, the mag- 
netic specific heat is expected to follow a C m (T) ~ 1/T 2 
form at temperatures large compared to the typical tempera- 
ture scale T p , the temperature at which the specific heat peaks, 
set by the interactions in these systems. This form was indeed 
verified in all our simulation results. In contrast, all the exper- 
imental C m (T) data decrease at T > 5 K significantly faster 
and are obviously not in agreement with this necessary 1/T 2 
high-temperature form. 

We believe this fast drop-off in experiment is likely due to 
the over-subtraction of the lattice contribution to the total spe- 
cific heat at these temperatures. The usual methods for carry- 
ing out such a subtraction rely on an estimated Debye contri- 
bution for the acoustic phonons. For example, by considering 
the temperature range of 10 K < T < 20 K, one might try to 
fit the total specific heat to the form C to t a i(T) = A/T 2 + BT 3 , 
where the 1/T 2 part comes from the aforementioned magnetic 
contribution and T 3 part is the Debye phonon contribution. 
Unfortunately, for T > 10 K, background contributions from 
other components of the experimental setup become signifi- 
cant. In particular, we note that in order to facilitate thermal 
conduction in the measurements, Ag powder was mixed into 
the spin ice powder. At these higher temperatures, the specific 
heat contribution from the Ag powder component becomes 
larger than the magnetic component that we are trying to iso- 
late. Fitting the phonon contribution with all these high tem- 
perature background contributions embeds errors in the A and 
B fitting parameters, which then causes an over-subtraction 
for the magnetic specific heat C m (T) at 5 K < T < 10 K. 



2. Low temperature regime 

We now turn to the low temperature regime of the C m (T) 
curves, below the prominent peak at T = T p , with T p ~ 1 
K for Dya-^Y^C^ and T p ~ 1.9 K for Hc^-zY^CV 
In particular, we discuss the minima found in the simulation 
results for all dilution levels (including x = 0, although in this 
case the minimum is more subtl o 28 ' 29 ) in both the Dy and Ho 
spin ices (see solid curves in insets in Fig. |2). As discussed in 
Subsection IIV Bl below, the integrated entropy of the system 
is highly dependent on the C m (T) results at low temperatures 

since dS = %^dT. 
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FIG. 2: (Color online) Comparison of the magnetic specific heat, C m (T), between Monte Carlo simulations and experiments. Black open 
circles are for Dy2-:c Y^T^C^ experiment, solid black curves are for Dy2-x Y^T^Ch simulations. Red open squares are for H02-H Y^Ti207 
experiment, and solid red curves are for Ho2-a: Y^TiaOy simulations. Insets show an enlargement around the Schottky peak at T p , arising from 
the formation of the spin ice state. The horizontal blue arrows indicate location of C m (T) minima that may be occuring in Ho2- x Y x Ti20r. 



It is known that in simulations of the undiluted dipolar spin 
ice modelj22i22 a C m (T) minimum arises from the develop- 
ment of extra correlations within the spin ice state caused by 
the dipolar interactions, with the system eventually undergo- 
ing a transition to long-range order at T c ~ 0.13D (T c ~ 0.18 
K, for the J±, D parameters appropriate for Fk^T^Oyi 28 ' 29 ) 
For such minima to be found in undiluted spin ice simula- 
tions, collective spin update algorithms (loop moves discussed 
in Section ITlIBI ) have to be included. On the other hand, it is 
very difficult for experiments to display such a C m (T) mini- 
mum and the long-range order transition, due to the freezing 
of spins below a temperature T ~ 0.5 KM 

For the diluted systems, the existence of the minima in our 
simulation suggests that a dynamical arrest similar to the one 
in the undiluted systems does occur. Indeed, as discussed 
in Section III B, equilibrium in the simulations cannot be 
achieved without using collective update algorithms, further 
supplemented by parallel tempering. For Dy2- a: Y 2 ;Ti207, 
having used a 3 He cryostat (See Section II B), the experiments 



stop at temperatures just above the simulation-predicted min- 
ima. For Ho2- 2; Y 2 ;Ti207, the C m (T) minima are perhaps ex- 
perimentally observed (see horizontal blue arrows in the in- 
sets of Fig. [2j, although the experimental data points below 
the minima do not agree very well with the simulation results. 
In this case, one should be warned that there is a large nuclear 
contribution at T < 0.5 K for I-^T^Ot 2 ^ Even though this 
nuclear component has been subtracted (see Section []]}, its 
existence nevertheless complicates the possible experimental 
observation of the minima in the magnetic-only part, C m (T), 
of the total specific heat C(T). 

While the present experimental data do not allow for a con- 
vincing observation of the minima in C m (T), we unquestion- 
ably find them in the Monte Carlo simulations of the micro- 
scopic DISMs. The minima observed in the simulations of the 
diluted DSIMs are significantly different from the ones in the 
undiluted variants i 28 ' 29 ! 43 Upon dilution, the C m (T) minimum 
acquires a significant value, as seen in Fig. [2] Furthermore, 
the broad specific heat peak at T p (x), which signals the devel- 
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opment of ice rule correlations as in the undiluted Dy and Ho 
spin ices, is less well defined in presence of dilution. For ex- 
ample, for x = 1.7, the peak is more that of a wiggly feature, 
on the rising C m (T) curve as T approaches zero, rather than 
a well-defined peak. Indeed, at such a high dilution, the ice 
rules are marginally enforced and the C m (T) peak associated 
with the development of ice rules fulfilling tetrahedra is not 
very prominent. As discussed further in Subsection IIV B I re- 
garding the determination of the residual entropy S Tes (To) at 
a low temperature T , the behavior of the C m (T) curves sug- 
gests that the residual entropy concept employed for undiluted 
spin ices cannot be readily discussed without a specification 
of the lowest temperature To at which (equilibrated) experi- 
mental data are obtained. 

To sum up, there exist significant systematic experimen- 
tal difficulties in determining the magnetic-only contribution 
to the specific heat, C m (T), in the high temperature regime 
(T > 5 K). For the low temperature regime (T < 0.5 K), 
in contrast to the undiluted case, the C m (T) curves from 
our simulations display clear minima with significant C m (T) 
values. On the experimental front, these minima may be 
marginally observed in Ho2_ x Y x Ti207 (x — 0.4,0.8,1.2), 
but are not observed in Dy2_ :c Y; ! ;Ti207. At the same time, 
the very good agreement between the experimental and Monte 
Carlo C m (T) for both materials (for x up to x = 1.8 for 
Dy 2 - x Y x Ti 2 7 ) and for 0.5 K < T < 5 K seemingly vin- 
dicates the applicability of a simple site-diluted version of the 
DSIM to describe Dya-^Y^Oy and ^-xY^TiaCV We 
thus take the following approach. Having demonstrated good 
agreement between experiments and models in the tempera- 
ture range T ~ [0.4 K - 5 K] for both Dy^YzTiaOy and 
Ho2- 2 ;Y x Ti207, in order to remedy the aforementioned ex- 
perimental caveats, we henceforth only consider the simula- 
tion data of Eq. (fJJ to expose accurately what is the x depen- 
dence of the low-temperature residual entropy, 5 res (To) of the 
Dy2_ a; Y a ;Ti207 and Ho2- a; Y :r Ti207 diluted dipolar spin ice 
materials. 



B. Non-monotonic Residual Entropy 

Since Eq. (f2]i is an Ising model, the entropy at infinite tem- 
perature per mole of spin is Rln 2. Thus the residual entropy 
at a given temperature Tq can be written as 



/>OC 

S res (T )=Rln2- / 

J T 



C m (T) 
T 



dT 



(4) 



We plot S rcs (T ) obtained from the Monte Carlo simulations 
for different choices of To, where the integration to infinite 
temperature are done by fitting the C m (T) curves at high tem- 
peratures (> 10 K) to the 1/T 2 form. 

The results from these Monte Carlo determinations of 
the residual entropy, 5 rcs (7o) are shown in Fig. [3] for both 
Dy2_ a; Y a ;Ti207 and Hc^-xYxTi^Oy. We confirm the pre- 
vious observation made by Ke et al. in Ref. l23ll that there 
does exist (i) a systematic non-monotonic x dependence of 
S res (T ) and (ii) that there is a difference in S rcs (T ) between 
the two materials. The main new result here is that, thanks to 
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FIG. 3: (Color online) Residual entropy determined from Monte 
Carlo simulations for both Dy2- a: Y a: Ti207 and fk^-iYzT^O- 
with different low temperature limits To. The dotted black curve 
shows SVcs given by the generalized Pauling's argument (gPa). 



the ability of the Monte Carlo simulations to provide accurate 
C m (T) data for T < 0.5 K and T > 10 K ranges, we can 
now robustly expose both the x dependence and the materi- 
als dependence of S rcs . Supplementing the previous report^ 
we are now also uncovering the importance of specifying the 
base temperature T used in the determination of S IOS (T ). 
Such a need to specify To does not arise in previous work 
on undiluted Dy2Ti2C<7 and Ho2Ti2C»7 because C m (T) prac- 
tically drops to zero near T ~ 0.4 K and S ICS remains close to 
the Pauling value for C m (T)/T integrated upward anywhere 
from 0.4 K ± 0.1 K. In particular, as a final and crucial ob- 
servation, we note that for all values of x and for a given To, 
S les (x) is Zower for Ho2- a; Y a ;Ti207 than for Dy 2-^3^2(1)7, 
in contrast to the conclusion that was reached in Ref. rf23ll and 
reproduced in Fig. Q] 

To reiterate, as can be seen in Fig. [3] the results of the resid- 
ual entropy for the diluted (x > 0) DSIM depend strongly on 
the choice of To, in contrast to the undiluted case {x = 0), 
in which the S rcs (To) for different T)S almost collapse onto 
the calculation of the Pauling's entropy, (R/2) ln(3/2). For 
x = 0, the collapse of the S rcs (To) for different To's is the 
manifestation of the projective equivalence^, which states 
that the quasi-ground state properties of the DSIM can be de- 
scribed by an effective nearest-neighbor spin ice model up to 
corrections falling off as 1/r 5 . But for x > 0, the To depen- 
dence suggests the failure of the projective equivalence upon 
dilution. 

The overall non-monotonic trend of the entropy from the 
generalized Pauling's argument being in rough qualitative 
agreement with the results for the real materials suggests 
a remnant of the diluted nearest-neighbor spin ice model 
physics in the diluted DSIMs. Yet, the two materials, be- 
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cause of their different magnetic interactions, display distinct 
5 res (x, Tq). Specifically, the two materials possess differ- 
ent energy scales for their dipolar interactions, D, relative 
to the nearest-neighbor energy scale, J\ (see Eq. (f2]i). Thus, 
the higher overall temperature scale for the formation of the 
spin ice state in Woi-x*L x^^O-; compared to Dy2- a jY 2; Ti207 
results in a residual entropy S rcs (To,x) for rk^-zY^^Ch 
lower than for Dy2_ a; Y 2 ;Ti207 for all x and for a given Tq. 
However, a choice of To that varies for different values of x 
for a given compound will lead to a less smooth 5 res (T , a;) 
evolution than the one seen in Fig.[3](see Fig.Q]). 



C. Large Level of Dilution 

It is perhaps remarkable that the nice agreement found be- 
tween Monte Carlo simulations and experiments shown in 
Fig. |2] for Dy2-a;Y x Ti207 for < x < 1.7 disappears 
abruptly and essentially completely going from x = 1.7 to 
x = 1.8 and x = 1.9 (see Fig. [2j. The only similarity left is 
that both Monte Carlo and experimental C m (T) data show a 
small low-temperature hump at a temperature T ~ 0.8 K that 
somewhat agrees between Monte Carlo and experiments (see 
insets of Fig.|2]for x = 1.8 and x = 1.9, which are further re- 
produced in Fig. |4j». This figure further illustrates that despite 
the large dilution of magnetic ions for x = 1.8 and x = 1.9, 
finite size effects remain negligible. We are thus rather confi- 
dent that the discrepancy between simulation and experimen- 
tal results does not arise from computational pitfalls, but is a 
genuine physical difference. 

Presently, we do not have a good suggestion as to what may 
cause such a sudden (in terms of "just" going from x = 1.7 
to x = 1.8) and large discrepancy between experiments and 
Monte Carlo data. A possible mechanism includes the devel- 
opment of a dipolar Ising spin glass stat o 56,57 inhibiting ther- 
mal equilibrium in the experiments, though that should not be 
at play at temperatures as high as 1 K. Another possibility in- 
cludes a significant random local lattice distortion developing 
upon reaching large levels of dilution. This would affect the 
J„ couplings and the crystal field, hence the magnetic moment 
fi and the coupling D compared to the values determined for 
x = 0. A third possibility is that of a highly uneven distri- 
bution of the magnetic ions as x — > 2. These last two pos- 
sibilities seem rather unlikely given the close ionic radius of 
Y 3+ with Dy 3+ and Ho 3+ and the solid solution that exist in 
the whole x G [0,2] range. More experiments are definitely 
required to understand the x — >• 2 behavior of diluted spin ice 
materials. 



V. CONCLUSION 

In this paper we have reported results from Monte Carlo 
simulations of a site-diluted version of the dipolar spin 
ice model (DSIM) given by Eq. [2] for Dy2_ :c Y. r Ti207 and 
Ho2-a;Y x Ti207. A close match between simulation results 
and experiments in the temperature range 0.5 K < T < 5J 
was found up to, and including, x = 1.7 (85% magnetic ions 
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FIG. 4: (Color online) Comparison of Monte Carlo specific heat with 
experimental results for Dy2- X Y^T^Ot for sizes L — 3,4,5 for 
x = 1.8 (top panel) and x = 1.9 (bottom panel). 



diluted) for Dy2-xYj;Ti207. This good agreement between 
simulations and experiments validates further the underlying 
dipolar spin ice models for these two compounds i^il 3 . 

The non-mono tonicity of the residual entropy as a func- 
tion of dilution levels, S ICS (T , x), is confirmed to originate 
from the material-specific spin-spin interactions themselves, 
namely the relative strength of the dipolar interactions with re- 
spect to the (mostly) nearest-neighbor exchange coupling J±. 
Furthermore, despite the importance of specifying the base- 
temperature To from which thermodynamic integration of the 
magnetic specific heat C m (T)/T is carried out, 5 res (To,a;) 
is nevertheless found to be roughly qualitatively described by 
the generalized Pauling's (gPa) estimate. In summary, the dif- 
ference in the residual entropy S ycs between Dy2- a: Y 2 ;Ti207 
and Ho2-a;Y x Ti207, as well as with the gPa, have been re- 
solved in the present work. 

Encouraged by the robustness of the site-diluted dipolar 
spin ice model to describe the experimental observations for 
temperatures higher than 0.5 K or so, we hope that our work 
will stimulate further experimental investigations and theoret- 
ical studies of spin ice materials at T < 0.5 K, in particular 
in the context of evincing a possible transition to long range 
or( j erj 28 i 29 j t wou j c j b e interesting to explore further the highly 
diluted regime of Dy2_ 2 :Y a ;Ti207 {x > 1.8) to clarify the ori- 
gin of the discrepancy between experimental and Monte Carlo 
specific heat data in that regime. It might also be interesting 
to explore the possibility of a dipolar Ising spin glass state in 
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the highly diluted regime of spin ice material s^Sai 
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